Method, device and computer program for monitoring a rotating machine of an aircraft

ABSTRACT

The invention relates to a method ( 1 ) for monitoring a rotating machine ( 100 ) of an aircraft, wherein a measurement signal is acquired from the rotating machine. According to the invention, instantaneous frequencies (f K (t)) of sinusoidal components of the measurement signal are estimated, and, using a computing module ( 12 ), a plurality of successive iterations are carried out in each of which: complex envelopes of the components are updated (C 1 ), parameters of a model of a noise of the signal are updated (C 21 ) on the basis of the envelopes, whether the model has converged from the preceding iteration to the present iteration is tested (C 4 ), with a view to: o if not, carrying out a new iteration, o if so, performing a computation (D) of the complex envelopes on the basis of the iterations that have been carried out.

The invention relates to a method, a device and a computer program for monitoring at least one rotating machine of an aircraft. In particular, it can be a rotating machine for propulsion of the aircraft.

The invention relates to the field of monitoring rotating machines by analysis of vibratory signals. Elements in rotation produce mechanical vibrations which can result from imbalance, eccentricity of mass, worn bearings, twisted shafts, misaligned shafts, or similar. The consequence of manifestations of imbalance is to create vibrations of cyclic character whereof the periodicity depends on the component and the phenomenon affected. For example, a defect in imbalance of a rotor generates a purely sinusoidal component the instantaneous frequency of which is equal to that of the rotation and amplitude proportional to the square of the speed of rotation. Another example is the defect in parallel alignment of a shaft consequently generating sinusoidal vibration at a frequency double that of the shaft. Such defects, and others, can cause reduction in the service life of the components of the rotating machine or untimely stopping of the mechanical system in the event of damage (for example if the shaft is broken). These defects can also affect other elements of the system. For example, a defect in alignment causes premature wear of bearings, joints, shafts and couplings in the kinematic chain.

For these reasons, it is often necessary to control these defects by monitoring the vibratory signal. The operator is often interested in evaluating the level of imbalance or poor alignment of the shaft during a run-up (from 0 speed to operating speed) and also in a drop on speed. It is also preferable to know if the rotor passes through one or more resonances from its start-up to its normal operating speed. A tracking filter can ensure this function by calculating the amplitude and phase of a sinusoidal component of interest during a run-up and a shut-down. A tracking filter is equivalent to a very narrowband filter whereof the central frequency follows the frequency in rotation of the component of interest. The latter is often proportional to the rotation frequency of the engine shaft.

For example, if the operator wants to control the alignment of a rotor of a rotating machine, a tracking filter with a central frequency equal to twice the rotation frequency of the shaft can be applied to extract the component of interest. The amplitude and the instantaneous phase of this component are valuable diagnostic indicators. The amplitude can indicate the level and the severity of the alignment defect (in some cases this can be tolerable) or the existence of resonance. On the other hand, a sudden change in phase may indicate passing through a resonance.

The first tracking filters were implemented by analogue electrical circuits, such as in U.S. Pat. Nos. 3,307,408A and 3,938,394A. The problems of analogue filters are the complexity of electrical circuits which often need considerable power consumption. For this, digital filters are preferred today. A tracking filter known and used widely in the field of rotating machines is the Vold-Kalman filter. The design of this filter is based on a physical model given to the vibratory signal. In fact, vibratory signals originating from rotating mechanical members can be modelled as a sum of harmonic components embedded in supplementary noise. Each component can be shown as a signal modulated in phase and amplitude by a low-frequency basis. Slow modulations in amplitude and phase are physically connected to the fact that a change in speed of the machine has the effect of exciting different regions of the transfer function of the system. This results in a change in gain and phase accompanied by the change in speed. The analytical expression of this known model is given by:

[Math.1] $\begin{matrix} {{y(n)} = {{\sum\limits_{k}{{a_{k}(n)}{\exp\left( {j{\phi_{k}(n)}} \right)}}} + {{noise}(n)}}} & (1) \end{matrix}$ [Math.2] $\left\{ \begin{matrix} {{y(n)}:{observed}{signal}} \\ {{\phi_{k}(n)}:{instantaneous}{phase}} \\ {{a_{k}(n)}:{complex}{envelope}} \end{matrix} \right.$

The instantaneous frequency f_(k)(n) describes the variation in frequential content over time. It is equal to the derivative of the instantaneous phase according to the following equation:

$\begin{matrix} {{f_{k}(n)} = {\frac{Fs}{2\pi}\left( {{\phi_{k}(n)} - {\phi_{k}\left( {n - 1} \right)}} \right)}} & \left\lbrack {{Math}.3} \right\rbrack \end{matrix}$

where Fs is the sampling frequency. It is important to note that the information of the amplitude and of the phase of the k-th component are respectively the amplitude and the phase of the associated envelope a_(k)(n).

The Vold-Kalman method proposes estimating the complex envelopes a_(k)(n) associated with the instantaneous frequencies f_(k)(n) from two equations. The first is called the equation of data stipulating that each component to be tracked is equal to the signal observed with a supplementary error δ_(k) (n) which must be as minimal as possible:

[Math. 4]

y(n)−a _(k)(n)exp(jϕ _(k)(n))=δ_(k)(n)

The Vold-Kalman method also needs information a priori on the component to be tracked, which will be integrated into the structural equation. The complex envelope to be estimated a_(k)(n) is a low-frequency signal modulating the carrier exp(jϕ_(k)(n)), implying a regular and smooth structure of the envelope. The basic idea behind the Vold-Kalman method is that of introducing local restrictions favouring this property. This leads to a second equation, called the structural equation, which profits from the discrete differences to act as a low pass filter:

[Math. 5]

∇^(p+1) a _(k)(n)=∈_(k)(n)

where ∇^(p) is the operator of difference in order p (p also designates the order of the filter) and ε_(k)(n) is a term of error which must be minimised to ensure the regularity of the envelope. In general, the orders 1 and 2 are the most used and the structural equations from which they derive are then given respectively by:

[Math. 6]

(p=1): a _(k)(n−1)−2a _(k)(n)+a _(k)(n+1)=∈_(k)(n)

[Math. 7]

(p=2): a _(k)(n−2)−3a _(k)(n−1)+3a _(k)(n)+a _(k)(n+1)=∈_(k)(n)

The selectivity of the filter increases for intersecting orders.

Given the structural equation and that of data, the estimator in terms of least squares of the complex envelope

[Math. 8]

a _(k)=[a _(k)(1), . . . ,a _(k)(N)]^(T)

is equivalent to minimising the norm

[Math. 9]

₂ ²

of the errors according to the following equation:

$\begin{matrix} {{\min\limits_{a_{k}}{\delta_{k}}^{2}} + {r_{k}{\epsilon_{k}}^{2}}} & \left\lbrack {{Math}.10} \right\rbrack \end{matrix}$ where $\begin{matrix} {{\delta_{k} = \left\lbrack {{\delta_{k}(1)},\ldots,{\delta_{k}(N)}} \right\rbrack^{T}},} & \left\lbrack {{Math}.11} \right\rbrack \end{matrix}$ $\begin{matrix} {\epsilon_{k} = \left\lbrack {{\epsilon_{k}(1)},\ldots,{\epsilon_{k}(N)}} \right\rbrack^{T}} & \left\lbrack {{Math}.12} \right\rbrack \end{matrix}$

and

r_(k) are positive weights to be fixed by the operator to control the degree of regularity of the envelope. The solution to this problem is explicit and is given in matrix form for each component k by the following equation:

[Math. 13]

=(I+D _(k) ^(T) D _(k))⁻¹ C _(k) ^(T) y

with

[Math. 14]

y=[y(1), . . . ,y(N)]T

[Math. 15]

D _(k)=√{square root over (r _(k))}A _(p),

where A_(p) is the matrix of discrete differences of order p and C_(k) is the diagonal matrix formed by the terms

[Math. 16]

exp(jϕ _(k)(1)), . . . , exp(jϕ _(k)(N)).

A first possibility for estimating several components

[Math. 17]

a ₁(n), . . . ,a _(K)(n)

would be to proceed sequentially and independently with estimation of each component, that is:

[Math. 18]

∀k∈{1, . . . K},

=(I+D _(k) ^(T) D _(k))⁻¹ C _(k) ^(T) y

The latter is known as mono-component estimation.

A second possibility would be to jointly estimate these components. In this second case, the equation of data is given as:

$\begin{matrix} {{{y(n)} - {\sum\limits_{k = 1}^{K}{{a_{k}(n)}{\exp\left( {j{\phi_{k}(n)}} \right)}}}} = {\delta(n)}} & \left\lbrack {{Math}.19} \right\rbrack \end{matrix}$

and the estimator in the sense of least squares of the vector of the complex envelopes

[Math. 20]

a=[a ₁ ^(T) . . . ,a _(K) ^(T)]^(T),

is the solution to the problem:

$\begin{matrix} {{\min\limits_{a}{\delta }^{2}} + {\sum\limits_{k = 1}^{K}{r_{k}{\epsilon_{k}}^{2}}}} & \left\lbrack {{Math}.21} \right\rbrack \end{matrix}$

which is given in matrix form by:

[Math. 22]

â=(C ^(T) C+D ^(T) D)⁻¹ C ^(T) y

with

[Math. 23]

=[C ₁ , . . . ,C _(K)]

and

D is the block-diagonal matrix formed by the sub-matrices D₁; . . . ; D_(k).

This estimation is multi-component.

However, neither of these two possibilities using the given Vold-Kalman method is satisfactory in the field of aeronautics.

There is a first problem not resolved by this known Vold-Kalman method, in the case of non-stationary and/or non-Gaussian noise. In fact, the Vold-Kalman method is optimal only assuming that the noise is Gaussian and stationary, which is not verified in aeronautics for two reasons:

-   -   Admittedly, aerodynamic noise can be approximated by Gaussian         noise but it is far from being stationary, especially in flight.     -   In the engine there are several rotating members which generate         periodic vibrations and can be modelled by sinusoidal models.         The components not considered in the structural equation will         act as noise adding to the aerodynamic noise and can distort         estimation of the envelopes of interest above all in the         presence of coupling between the components.

For these reasons, noise (which is what remains in the measured signal if the components of interest are subtracted), is composed of random noise of Gaussian type (stationary or not) and sinusoidal noise (residual interfering sinusoidal components which are not included in the model among the K components to be estimated). This sinusoidal noise is far from being Gaussian. Also, at variable speed (when the engine is rotating at a non-constant speed/charge), this sinusoidal noise is non-stationary.

To approach the assumption of Gaussian noise and attenuate artefacts at coupling instants, joint estimation (multi-component) of the signal of interest with all the other harmonics which intersect it or which are very close over a period of time is preferable. In fact, when the different interfering components are estimated jointly, the equation of data ensures that the total energy of the signal will be distributed between these components, and with the structural equation favouring the regularity of solutions, joint estimation best decouples the components the frequencies of which intersect.

Therefore, this first non-resolved problem creates several disadvantages:

-   -   The major disadvantage of multi-component estimation in the         Vold-Kalman method is a higher calculation relative to         mono-component estimation: the dimension of the matrix to be         inverted is K times greater. Taking all vibrations into         consideration, K is going to be substantial.     -   Most often, the operator is not interested in the results for         estimating interfering components. These components are         appropriately estimated to diminish the estimation artefacts of         the components of interest.     -   The total number of sinusoidal components in the signal measured         by the sensor is most often unknown.     -   To detect the presence of interfering components, the operator         must perform an inspection of the spectrogram of a measured         multi-component signal. Also, it can happen that the         instantaneous frequencies of these interfering harmonics are         unknown, requiring extra work in estimating these frequencies         prior to proceeding with estimating by Vold-Kalman, which can be         prohibitive in the event where there is a large number of         interfering components.     -   To attenuate these artefacts, a large value of r_(k). can be         used, but this supposes that the coupling moments a priori are         known, which is not the case in practice. Placing a high value         of r_(k) on the whole signal is not a reliable solution, as         explained hereinbelow in reference to FIG. 4 d , as that leads         to underestimation of the amplitude of the envelope.

There is a second problem not resolved by this known Vold-Kalman method, relating to the difficulty of adjusting the parameter r_(k). Known methods need manual adjusting of parameters r_(k), or, as an equivalent, of the bandwidth associated with each sinusoid. The choice of these parameters r_(k) significantly affects the results of estimation and therefore a good choice is required. The quality of the estimation is closely tied to the choice of parameter r_(k). The parameter r_(k) can be interpreted statistically as the ratio of variances of the two centred random Gaussian vectors δ_(k) and ε_(k) respectively. It follows that the greater r_(k) is, the smaller the variance of ε_(k) and therefore the solution is increasingly influenced by the structural equation resulting from a very smooth envelope. This also means that the greater r_(k) is, the narrower the bandwidth of the filter associated with the k-th frequency and therefore the greater the selectivity. The optimal value of r_(k) depends on the signal level over noise. The greater the noise, the greater the need to use a large value of r_(k) so the noise can be suppressed. However, adjusting this parameter r_(k) remains difficult for several reasons:

-   -   If noise is non-stationary, as is the case of a rotating         aircraft machine, especially for propulsion of the aircraft (the         level of this noise can change in flight, for example), use of         the same value of r_(k) is not optimal.     -   The use of a very big value of r_(k) can cause digital problems.         In fact, the matricial product D_(k) ^(T)D_(k) gives a         semi-defined positive matrix of which the maximal eigen value         increases linearly with r_(k). In mono-component estimation, the         matrix I+D_(k) ^(T)D_(k) is a positive defined symmetrical         matrix (and therefore invertible) the minimal eigen value of         which is equal to 1. In multi-component estimation, C^(T)C is a         positive semi-defined matrix, so the minimal eigen value of         C^(T)C+D^(T)D can be very low. Therefore the values of the         weights r_(k) must not exceed a certain limit, beyond which the         matrix I+D_(k) ^(T)D_(k) (OU C^(T)C+D^(T)D) in the case of         multi-component estimation) can become very badly conditioned,         resulting in a problem of inversion and amplification of digital         errors. FIG. 1 shows the evolution of the signal-to-noise ratio         SNR of the component estimated in a noisy synthetic         mono-component signal as a function of the selected value of the         parameter r=r_(k). It can be seen that beyond a certain limit of         the parameter r_(k) dependent on the order of the filter, this         signal-to-noise ratio SNR of the estimated component drops         sharply due to poor conditioning of the matrix to be inverted         and that consequently the performance of the method is degraded.     -   The problem of adjusting the parameter r_(k) is more difficult         in the event of many harmonics (that is, K is large) are to be         tracked. In fact, the value of r_(k) does not depend only on the         level of the noise but also on the statistics of the components         to be estimated. In this way considering the same value for all         the r_(k) is generally not a reliable solution, given that the         different components to be estimated can have different energy         ranges and different variation rates. This problem can be         encountered for example if the aim is to estimate components         linked to the meshing of rotating cogged wheels in the rotating         machine and the modulations. In fact, the modulations generally         have low energy relative to the component linked to meshing.     -   Several defects in the rotating machine have an effect of         generating sinusoidal components at known frequencies. The         problem amounts to tracking this frequency and analysing its         amplitude. In the event where there is no defect, the amplitude         should be low, but poor adjusting of the Vold-Kalman method can         produce a dummy component whereof the amplitude is equal to the         smoothed background noise, which can cause errors in diagnosis.         Comparison of the amplitude estimated with the background noise         can resolve the presence of these false positives. But this         presupposes that the noise is known.

Given these observations, it appears that the problem of adjusting the parameters r_(k) is not obvious to resolve and can affect the quality of the filtering, which consequently compromises the efficacy of the Vold-Kalman method.

The invention aims to provide a method, a device and a computer program for monitoring at least one rotating machine of an aircraft, which resolves both the abovementioned problems for coupling components and adjusting the parameter of the Vold-Kalman method.

For this purpose, a first subject matter of the invention is a monitoring method of at least one rotating machine of an aircraft, a method in which

-   -   at least one measurement signal of the rotating machine is         acquired by an acquisition module,

characterised in that

-   -   instantaneous frequencies of sinusoidal components of the         measurement signal are estimated by an estimation module,     -   several successive iterations are executed by a calculation         module, in each of which:         -   complex envelopes of the sinusoidal components are updated,         -   parameters of a noise model of the measurement signal are             updated from the updated complex envelopes,         -   a convergence test is performed as to whether the model             converges from the preceding iteration to the current             iteration, to:             -   in the negative, redo a new iteration,             -   in the affirmative, perform a calculation of the complex                 envelopes from the iterations carried out.

The invention offers the possibility especially of resolving the second abovementioned problem of adjusting the parameter of the Vold-Kalman method. The invention allows automating of the estimation problem, and generalising to non-stationary interferences in particular the coupling of components and the non-stationarity of aerodynamic noise.

The invention therefore provides a method and a device for extraction or monitoring of one or more sinusoidal vibratory components having a variable instantaneous speed. By comparison with the Vold-Kalman method, the invention differs by way of the following advantages:

-   -   The problem is formulated probabilistically. The outlets of the         method are no longer specific estimators but laws describing         possible solutions, which makes it possible to have not only         specific estimators but also quantification of the confidence in         the estimation.     -   The invention allows automation of the estimation by calculating         some r_(k). jointly with the envelopes of the estimators. This         avoids adjusting manually or heuristically of the parameter by         the operators, and maximises the efficacy of the filtering.

According to an embodiment of the invention, all the unknown variables of the problem: the complex envelopes, the parameters of the noise model, the regularisation parameters, are estimated from their conditional laws a posteriori given the measurements. The parameters of these laws a posteriori are calculated by the calculation module at each iteration. As is known to the skilled person in general, a law a priori is defined by the fact that it is selected by the user so as to incorporate the preferred characteristics of the planned solution. As is known to the skilled person in general, a law a posteriori is defined by the fact that it takes observations into account.

According to an embodiment of the invention, said calculation of the complex envelopes in said affirmative is performed from the average of the law estimated.

According to an embodiment of the invention, the noise model of the measurement signal comprises a structured interfering noise, whereof the distribution of probability a priori of its spectrum has a heavy-tailed law calculated by the calculation module at each iteration. Therefore, according to an embodiment of the invention, the noise model of the measurement signal comprises a structured interfering noise, whereof the distribution of probability of its Fourier transform belongs to the class of heavier-tailed distributions than the normal law. This favours parsimony of the structured sinusoidal noise in the Fourier field.

According to an embodiment, the heavy-tailed law of the spectrum of the noise measurement is of Bernoulli-Gaussian type. According to an embodiment, the heavy-tailed law is given by a barycentre between a first distribution of Gaussian type calculated by the calculation module at each iteration and weighted by a weight, β or β_(j), non-zero and calculated by the calculation module at each iteration, and a second Dirac distribution calculated by the calculation module at each iteration and weighted by a weight 1−β or 1−β_(j), non-zero and calculated by the calculation module at each iteration.

According to an embodiment of the invention, the noise model of the measurement signal comprises a non-structured interfering noise Gaussian whereof the parameters are calculated by the calculation module at each iteration. Therefore, according to an embodiment of the invention the noise model of the measurement signal comprises non-structured interfering noise whereof the distribution of probability a priori is given by a third distribution of Gaussian type calculated by the calculation module at each iteration.

According to an embodiment of the invention, the noise model of the measurement signal resulting from a structured interfering noise model and from a non-structured interfering noise model is the convolution of a Gaussian law and of a heavy-tailed law whereof the parameters are calculated by the calculation module at each iteration. According to an embodiment, the heavy-tailed law of the spectrum of the noise measurement is of Bernoulli-Gaussian type. According to an embodiment, the convolution is a third distribution of Gaussian type calculated by the calculation module at each iteration and weighted by a non-zero weight, 1−β or 1−β_(j), calculated by the calculation module at each iteration and a fourth distribution of Gaussian type calculated by the calculation module at each iteration and weighted by a non-zero weight, β or β_(j), of the presence of structured noise and calculated by the calculation module at each iteration, the third distribution and the fourth distribution having variances different to each other and calculated by the calculation module at each iteration.

According to an embodiment of the invention, the parameters of the noise model of the measurement signal comprise at least the weight β or β_(j) and/or 1−β or 1−β_(j).

According to an embodiment of the invention, the parameters of laws of structured noise and of non-structured noise are considered as random variables to which laws of probability are assigned conjugated with the noise model. According to an embodiment of the invention, a conjugated law of probability a priori is allocated to hyperparameters of laws a priori of the parameters of the noise model. According to an embodiment of the invention, the precision parameter of the Gaussian law of non-structured noise is modelled by a Gamma law. According to an embodiment of the invention, in the event of a Bernoulli-Gaussian law for structured noise the probability β of the presence of the structured noise is modelled by a uniform law a priori and the precision parameter of the Gaussian law of structured noise is modelled by a Gamma law a priori.

According to an embodiment of the invention, a Beta law of probability a posteriori is followed by the weight β or β_(j) and/or 1−β or 1−β_(j) and is calculated by the calculation module at each iteration.

According to an embodiment of the invention, a law of Gamma probability a priori is allocated to at least one τ_(j) and a law of Gamma probability a posteriori is calculated by the calculation module at each iteration, where 2τ_(j) ⁻¹ is the variance of the third distribution.

According to an embodiment of the invention, a law of Gamma probability a priori is allocated to at least one η_(j) and a law of Gamma probability a posteriori is calculated by the calculation module at each iteration, where

[Math. 24]

η_(j)=1/(ζ_(j) ⁻¹+τ_(j) ⁻¹)

is the variance of the fourth distribution,

2τ⁻¹ is the variance of the third distribution,

2ζ_(j) ⁻¹ is the variance of the first distribution of Gaussian type.

According to an embodiment of the invention, the complex envelopes to be estimated are modelled by a law ensuring regularity of the envelopes and depends on a regularisation parameter γ_(k) which plays the same role as the parameter r_(k) in the Vold-Kalman method.

$\begin{matrix} {{p\left( {a^{R},{a^{I}❘\gamma_{1}},\ldots,\gamma_{K}} \right)} = {{\prod\limits_{k = 1}^{K}{{p\left( {a_{k}^{R}❘\gamma_{k}} \right)}{p\left( {a_{k}^{I}❘\gamma_{k}} \right)}}} = {C_{a}{\prod\limits_{k = 1}^{K}{\gamma_{k}^{N}{\exp\left( {{- \gamma_{k}}{{A_{p}a_{k}}}^{2}} \right)}}}}}} & \left\lbrack {{Math}.25} \right\rbrack \end{matrix}$

where

[Math. 26]

C _(a)

is an constant independent from

[Math. 27]

γ₁, . . . ,γ_(K)

and

[Math. 28]

a

According to an embodiment of the invention, the parameter γ_(k) is considered as a variable random and modelled by a conjugated Gamma law.

According to an embodiment of the invention, an auxiliary variable is added to the model to decouple the different components to be estimated. The interest of this auxiliary variable is not physical but is for implementation reasons only.

According to an embodiment of the invention, the auxiliary decoupling variable is defined by a model having a distribution of probability of normal circularly symmetrical complex law type, which is calculated by the calculation module at each iteration.

According to an embodiment of the invention, the first parameters of the noise model of the measurement signal comprise at least one variance of the or of at least one of the distributions of Gaussian type.

According to an embodiment of the invention, the calculation module is configured, at each of the successive iterations carried out, to:

-   -   update the law of the complex envelopes from the parameters of         the noise model, the auxiliary variable, the regularisation         parameters and the observation,     -   update the parameters of the noise model from the updated         envelopes.

According to an embodiment of the invention, at each iteration regularisation parameters of the model are updated by the calculation module from the updated complex envelopes.

According to an embodiment of the invention, a Gamma law of probability a priori is allocated to at least one of the regularisation parameters and results from a Gamma law a posteriori whereof the parameters are calculated by the calculation module at each iteration.

According to an embodiment of the invention, at each iteration an auxiliary decoupling variable of the complex envelopes of the components between them in the model is updated by the calculation module, from the updated complex envelopes and from the parameters of the updated noise model.

According to an embodiment of the invention, the auxiliary decoupling variable is defined by a model having a distribution of probability of complex circularly symmetrical normal law type, which is calculated by the calculation module at each iteration.

According to an embodiment of the invention, the convergence test comprises testing by the calculation module at each current iteration that the difference between the absolute value of the envelope calculated from the current iteration and the absolute value of the envelope calculated from the preceding iteration is less than a prescribed threshold.

According to an embodiment of the invention, the complex envelopes having been calculated are analysed by the calculation module, and the communication of an alarm is triggered by the calculation module on at least one physical outlet as a function of the complex envelopes having been analysed.

A second subject matter of the invention is a monitoring device of at least one rotating machine of an aircraft for executing the method such as described hereinabove, the device comprising

-   -   a module for acquiring at least one measurement signal of the         rotating machine,

characterised in that the device comprises

-   -   a module for estimating instantaneous frequencies of sinusoidal         components of the measurement signal,     -   a calculation module configured, at each of several successive         iterations carried out, to:         -   update complex envelopes of the sinusoidal components,         -   update parameters of a noise model of the measurement signal             from the updated complex envelopes,         -   test whether the model converges from the preceding             iteration to the current iteration, to:             -   in the negative, redo a new iteration,             -   in the affirmative, perform a calculation of the complex                 envelopes from the iterations carried out.

According to an embodiment of the invention, the complex envelopes of the sinusoidal components are updated independently.

According to an embodiment of the invention, the calculation module is configured, at each of the successive iterations carried out, to:

-   -   update regularisation parameters of the model from the updated         complex envelopes

According to an embodiment of the invention, the calculation module is configured, at each of the successive iterations carried out, to:

-   -   update an auxiliary decoupling variable of the complex envelopes         of the components between them in the model, from the updated         complex envelopes and from the parameters of the updated noise         model.

According to an embodiment of the invention, the calculation module is configured to analyse the complex envelopes having been calculated, and is able to trigger the communication of an alarm on at least one physical outlet of the device as a function of the complex envelopes having been analysed.

A second subject matter of the invention is a computer program comprising code instructions for executing the monitoring method such as described hereinabove, when it is executed on a computer.

A third subject matter of the invention is an aircraft comprising at least one rotating machine and a monitoring device of the at least one rotating machine such as described hereinabove.

The invention will be better understood from the following description given solely by way of non-limiting example in reference to the figures of the appended drawings, in which:

FIG. 1 schematically illustrates the evolution of the signal-to-noise ratio of a component estimated according to a method of the prior art,

FIG. 2 illustrates an explanatory diagram of the device and of the monitoring method according to an embodiment of the invention,

FIG. 3 illustrates an explanatory flowchart of a part of the monitoring method according to an embodiment of the invention,

FIG. 4 a shows the initial spectrogram of a measured signal,

FIG. 4 b shows an example of an envelope having been calculated by the device and the monitoring method according to an embodiment of the invention,

FIG. 4 c shows the spectrogram of a residual signal obtained from FIGS. 4 a and 4 b by the device and the monitoring method according to an embodiment of the invention,

FIG. 4 d shows the spectrogram of a residual signal obtained by adjusting a method of the prior art,

FIG. 4 e shows the spectrogram of a residual signal obtained by other adjustment of a method of the prior art,

FIG. 4 f shows the spectrogram of a residual signal obtained by even other adjustment of a method of the prior art,

FIG. 5 illustrates a modular block diagram of the monitoring device according to an embodiment of the invention.

The monitoring device 1 of one or more rotating machines 100 embedded on an aircraft is described hereinbelow in reference to the figures. The rotating machine 100 can be a rotating propulsion machine of the aircraft, which can be one or more turbomachines, such as for example one or more turbojets or one or more turboprop. The aircraft can be a plane or a helicopter.

The monitoring device 1 performs the steps of the monitoring method of the rotating machine 100, described hereinbelow.

In FIGS. 2 and 5 , the monitoring device 1 comprises an acquisition module 10 for acquiring one or more measurement signals y(t) for measuring of one or more magnitudes on the rotating machine 100 while this rotating machine 100 is operating as a function of the time t. The measurement signal y(t) can be a measurement signal y(t) of vibrations or acoustics, for example of a rotor and/or a stator, the signal being due to rotation of the rotor relative to the stator. The measurement signal y(t) of vibrations can be measured on a helicopter on a gas generator and/or on a gearbox and/or on a meshing and/or on a bearing. The measurement signal y(t) of vibrations can be measured on a plane, on a stator of a turbomachine. The acquisition module 10 can comprise one or more sensors, embedded on or near the rotating machine 100 on the aircraft, for measuring the signal or signals y (t), which can comprise for example one or more measurement accelerometers of the vibratory signal y(t) and/or one or more measurement microphones of the vibratory signal y(t) acoustic, and/or one or more sensor (encoder or other) of rotation speed of a shaft or of a rotor of the rotating machine 100, or other.

According to an embodiment, the acquisition module 10 enables to measure and acquire a vibratory or acoustic signal generated by the rotating machine 100 operating. The sensor enables to obtain an analogue signal representative of the vibratory or acoustic signal which is connected to an acquisition chain 13 configured to provide the digital measurement signal y(t)=y (nTe) with T_(e) the non-zero sampling period, where T=nT_(e) represents different sampling instants, where n can be for example natural integer numbers. As a consequence, hereinbelow the variable n is a temporal variable. This acquisition chain 13 in this case comprises a conditioner, an anti-aliasing analogue filter, a sampler blocker and an analogue-digital converter.

The measurement signal y(t) comprises K sinusoidal components of interest x_(k)=a_(k)(n)exp(jΦ_(k)(n)), the sum of which is equal to x(n) according to the equation hereinbelow:

$\begin{matrix} {{x(n)} = {\sum\limits_{k = 1}^{K}{{a_{k}(n)}{\exp\left( {j{\phi_{k}(n)}} \right)}}}} & \left\lbrack {{Math}.29} \right\rbrack \end{matrix}$

I={1, . . . , K} designates the set of indices of components of interest of the measurement signal y(t).

Φ_(k)(n) designates the phase of each sinusoidal component of interest x_(k)=a_(k)(n)exp(jΦ_(k)(n)).

a_(k) (n) designates the complex amplitude (or complex envelope mentioned hereinbelow) of each component of interest a_(k)(n)exp(jΦ_(k)(n)) or the temporal derivative of any order of this component of interest a_(k) (n)exp(jΦ_(k) (n)).

In FIGS. 2 and 5 , the monitoring device 1 comprises an estimation module 11 for estimating several instantaneous frequencies f_(k)(t) of the sinusoidal components x_(k)=a_(k)(n)exp(jΦ_(k)(n)) of the measurement signal y(t), where k is a natural integer number varying from 1 to K, and where K is a prescribed natural integer number greater than 1. These instantaneous frequencies f_(k)(t) are frequencies to be monitored in the rotating machine 100 (for example the first order of the gas generator, of the free turbine, of the shaft of the reducer, the meshing frequencies of the gearbox and of the meshings in the starter, or other).

According to an embodiment of the invention, a rotation frequency f_(ref)(t) of a reference shaft of the rotating machine 100 (which can be for example the shaft of its rotor) can have been prescribed in advance in the estimation module 11. The rotation frequency f_(ref)(t) can be calculated from the speed signal measured by a rotation speed sensor forming part of the acquisition module 10 or from the vibratory or acoustic signal y(t) provided by the acquisition module 10. The estimation module 11 can calculate the instantaneous frequency f_(k)(t) of each sinusoidal component x_(k)=a_(k)(n)exp(jΦ_(k)(n)) as being proportional to the rotation frequency f_(ref)(t) according to the equation f_(k)(t)=c_(k)f_(ref)(t), where each proportionality coefficient c_(k) of the instantaneous frequency f_(k)(t) is real. Since the connection between the different components of the kinematic chain of the rotating machine 100 is rigid, the proportionality coefficients c_(k) can have been prescribed in advance in the estimation module 11. In FIG. 2 , the elements of the machine 100 which are monitored are specified in advance and their kinematic is prescribed in the estimation module 11, without any need for the user to give the frequencies f_(k)(t) in entries.

For example hereinbelow n is a natural integer number ranging from 1 to N, where N is a prescribed natural integer number greater than 1.

It is assumed that the measurement signal y(n) is equal to the sum of the K components of interest a_(k)(n)exp(jΦ_(k)(n)) to which the non-structured interfering signal (or noise) l(n) of Gaussian type (also called wideband) is added, and the structured interfering signal (or noise) e(n) (also called narrowband), according to the following equation:

y(n)=x(n)+e(n)+l(n)

According to an embodiment of the invention, the structured interfering signal (or noise) e(n) is defined as being the part of the interfering noise whereof the Fourier transform is discrete, that is this noise e(n) can be written as the sum of sinusoidal signals. In practice, this noise e(n) is constituted by the residual sinusoidal components not included in the signal x(n). The Fourier transform of the structured interfering noise e(n) is constituted by peaks each linked to a sinusoid. It is said that the spectrum of the structured noise e(n) is characterised by a parsimonious look.

According to an embodiment of the invention, the non-structured interfering signal (or noise) l(n) is defined as being that part of the interfering noise whereof the Fourier transform is not discrete. The non-structured interfering noise l(n) is constituted by the non-sinusoidal random components whereof the spectrum is not discrete. This non-structured interfering noise l(n) can be for example aerodynamic noise and is modelled by a Gaussian law.

According to the invention, e(n) incorporates a random signal following a law other than a normal law and must favour a parsimonious spectrum.

As a reminder and by definition in general, a normal or complex Gaussian law f(x)=NC(μ,σ²) in a real dimension is defined as the distribution of average μ and standard deviation a or variance σ² having the density of probability f(x) according to the following equation:

$\begin{matrix} {{f(x)} = {\frac{1}{\sigma\pi}e^{- {(\frac{x - \mu}{\sigma})}^{2}}}} & \left\lbrack {{Math}.30} \right\rbrack \end{matrix}$

The signal (or structured interfering noise) e(n) represents the residual harmonic signal, and can be noted in the following form:

$\begin{matrix} {{e(n)} = {\sum\limits_{k \in \overset{¯}{I}}{{a_{k}^{\prime}(n)}{\exp\left( {j{\phi_{k}(n)}} \right)}}}} & \left\lbrack {{Math}.31} \right\rbrack \end{matrix}$

where

[Math. 32]

Ī

is the set of indices of the sinusoidal components of the structured interfering signal e (n), which are each distinct from the indices of the set I.

The total interfering signal or total noise of the measurement signal y (t) is equal to the sum b (n) of the non-structured Gaussian interfering signal (or noise) l(n) and of the structured interfering signal (or noise) e(n) according to the following equation:

b(n)=e(n)+l(n)

Calculation of the sinusoidal components of interest a_(k)(n)exp(jΦ_(k)(n)) is of Bayesian type. Therefore, n, x(n), e(n) and l(n) are the temporal successive realisations of random variables the laws of which will be defined hereinbelow. Also, x (n), e(n) and l(n) are independent of each other.

Hereinbelow, the exponent R designates the real part and the exponent I designates the imaginary part (distinct from the above set I).

The vector y^(R) is defined as being the real observation vector, formed from the real N values y (n)^(R) of y (n) for the natural integer number n ranging from 1 to N, according to the following equation:

[Math. 33]

y ^(R)[y(1)^(R) , . . . ,y(n)^(R) , . . . ,y(N)^(R)]^(T)

where the exponent T designates the transpose.

The vector y^(I) is defined as being the imaginary vector of observation, formed by the N imaginary values y(n)^(I) of y (n) for n ranging from 1 to N, according to the following equation:

[Math. 34]

y ^(I)=[y(1)^(I) , . . . ,y(n)^(I) , . . . ,y(N)^(I)]^(T).

The vectors x^(R),x^(I), x_(k) ^(R), x_(k) ^(I), I^(R), I^(I),e^(R),e^(I) are also formed by the real or imaginary N values of x (n) or l(n) and e (n), for n ranging from 1 to N.

Noise Model

A law of observation of noise b(n) or noise model b(n) is defined hereinbelow. This noise model b(n) is calculated by a calculation module 12 (or monitoring unit 12) in FIG. 2 .

The calculation module 12 is digital (as opposed to analogue circuits). The calculation module 12 is automatic in the sense where it automatically performs the steps described hereinbelow. As shown in FIG. 5 , the calculation module 12 can be or comprise one or more calculators, and/or one or more computers, and/or one or more calculating machines, and/or one or more processors and/or one or more microprocessors. The calculation module 12 and/or the acquisition module 10 can comprise a computer program loaded in their memory to automatically execute the steps described hereinbelow. Also, the monitoring device 1 can comprise one or more physical outlets SP for presentation of information to present to a user the complex envelopes a_(k)(n) having been calculated by the module 12, wherein this physical outlet or these physical outlets SP can comprise a display screen for displaying these complex envelopes a_(k)(n) having been calculated (or their derivative of any order (amplitude, phase, amplitude as a function of speed, the estimated component or the estimated components, etc.)). The memory is provided for storing the outlets a_(k)(n) having been calculated by the calculation module 12.

Structured Noise Model e(n) The structured interfering noise model e(n) can be realised by one of the embodiments described hereinbelow.

The noise model e(n) is defined as being a structured interfering signal as follows. According to an embodiment of the invention, the spectrum of the structured interfering noise e(n) is defined by a model having a parsimonious spectrum, that is, a distribution concentrated at zero for favouring zero values with a large probability and having heavier tails than the normal law to authorise the presence of high values. For example, a good distribution of probability of the spectrum of structured interfering noise e(n) can be defined by a model being a mix of a finite/infinite number of different weighted distributions.

According to an embodiment of the invention, the spectrum of narrowband interfering noise is defined by a weighted sum of a distribution of probability of Gaussian type and of another distribution of probability of non-Gaussian type.

Of the distributions of infinite mixes, there is for example the Laplace law and the Student law currently used for modelling parsimonious signals and which can be expressed as infinite mixes of Gaussians with parameters of mixes following gamma laws and inverse gamma laws respectively. Of the distributions of finite mixes, there is for example the mix of two distributions: of a Bernoulli distribution and of a Gaussian distribution known as Bernoulli-Gaussian. In the case of the Bernouilli model, the spectrum of the structured interfering noise e(n) is defined by a model having a distribution of probability given by the barycentre between a distribution of Gaussian type NC(0, 2ζ⁻¹) weighted by a weight (probability) β of the presence of interfering noise e(n) and a Dirac distribution δ₀ weighted by a weight 1−β, for example as per the following equation:

[Math. 35]

∀n∈{1, . . . ,N}[F _(N) e]_(n)|β,ζ˜(1−β)δ₀ +βNC(0,2ζ⁻¹)

where

F_(N) designates the discrete Fourier transform operator of size N, this operator being divided by

[Math. 36]

√{square root over (N)},

δ₀ is the Dirac distribution positioned at the zero abscissa,

0≤β≤1, is the probability of the presence of the structured noise,

NC(0, 2ζ⁻¹) designates a normal complex distribution of zero average, whereof the variance is equal to 2ζ⁻¹. Therefore, the real part of this normal complex distribution NC(0, 2ζ⁻¹) is a Gaussian distribution of zero average and variance ζ⁻¹. The imaginary part of this normal complex distribution NC(0, 2ζ⁻¹) is a Gaussian distribution of zero average and variance ζ⁻¹.

According to an embodiment, β is different from zero.

Hereinabove, [F_(N)e]_(n) designates the discrete Fourier transform operator of size N, which is divided by

[Math. 37]

√{square root over (N)}

and is applied to e (n). Therefore, there is

$\begin{matrix} {F_{N} = {\frac{1}{\sqrt{N}}{DFT}_{N}}} & \left\lbrack {{Math}.38} \right\rbrack \end{matrix}$ And $\begin{matrix} {{F_{N}^{\top}F_{N}} = {{F_{N}F_{N}^{\top}} = {Id}}} & \left\lbrack {{Math}.39} \right\rbrack \end{matrix}$

where Id designates the matrix identity.

Unlike the known Vold-Kalman method, where the noise is assumed as a Gaussian centre, this embodiment easily takes structured noise into account in the model to improve the efficacy of the estimation at the time of coupling without need to estimate the interfering components.

In variable mode, the spectrum of e(n) can have wider peaks. Fourier transforms S_(j), limited to a temporal segment of size N_(j) for e(n) are defined as follows:

[Math. 40]

∀j∈{1, . . . ,J} ∀n∈{1, . . . ,N _(j)}[S _(j) e]_(n)|β_(j),ζ_(j)˜(1−β_(j))δ₀+β_(j) NC(0,2ζ_(j) ⁻¹)

with

[Math. 41]

S _(j) =F _(N) _(j) P _(j)

where P_(j) is a selection matrix of size

[Math. 42]

N _(j) ×N,

N_(j) is a prescribed natural integer number less than N (N_(j) is the length of the window),

$\begin{matrix} {F_{N_{j}} = {\frac{1}{\sqrt{N_{j}}}{DFT}_{N_{j}}}} & \left\lbrack {{Math}.43} \right\rbrack \end{matrix}$ with $\begin{matrix} {DFT}_{N_{j}} & \left\lbrack {{Math}.44} \right\rbrack \end{matrix}$

the discrete Fourier transform of size N_(j),

n is a natural integer number ranging from 1 to N_(j),

J is a prescribed natural integer number (number of temporal windows j), j is a natural integer number ranging from 1 to J,

NC(0, 2ζ⁻¹) designates a normal complex distribution of zero average, whereof the variance is equal to 2ζ⁻¹,

β_(j) is the weight of the presence of structured noise, allocated to NC(0, 2ζ⁻¹),

0≤β_(j)≤1,

1−β_(j) is the weight allocated to the Dirac distribution positioned at the zero abscissa. Fourier transforms S_(j) are therefore of shorter duration than F_(N).

If, also, the matrices P_(j) have disjointed supports, this gives

[Math. 45]

P _(j) P _(j) ^(T) =Id

and therefore

[Math. 46]

S _(j) S _(j) ^(T) =Id

The matrices P_(j) have disjointed supports if the sum of all the N, is equal to N:

$\begin{matrix} {{\sum\limits_{j = 1}^{J}N_{j}} = N} & \left\lbrack {{Math}.47} \right\rbrack \end{matrix}$

The selection matrix P_(j) is such that it selects N_(j) values belonging to the window j of a vector of size N.

P_(j) is a matrix containing N_(j) lines and N columns, it is defined by:

For j=1:

[Math. 48]

∀i ₀∈{1, . . . ,N _(j) }i ₁∈{1, . . . ,N}P _(j)(i ₀ ,i ₁)=1 if i ₁ =i ₀ if not P _(j)(i ₀ ,i ₁)=0.

For j>1

[Math. 49]

∀i ₀∈{1, . . . ,N _(j) }i ₁∈{1, . . . ,N}P _(j)∈(i ₀ ,i ₁)=1 if i ₁ =i ₀+Σ_(j) ₀ ₌₀ ^(j−1) N _(j) ₀ if not P _(j)(i ₀ ,i ₁)=0.

Unlike the known Vold-Kalman method where noise is assumed to be stationary, which is rarely the case of aerodynamic noise, this embodiment takes into account the non-stationary aspect of noise, and does this by analysing the signal via short-term Fourier transforms S_(j) mentioned hereinabove. The time axis is divided into J intervals or windows. The parameters of the model can change from one window to the other. Given that the parameters of these laws are unknown (the level of aerodynamic noise is unknown and the coupling instants between the components and the amplitude of the interfering components are unknown), they are also going to be adapted automatically. Laws a priori are selected to model information a priori on these parameters. In general, if there is no information a priori laws called non-informative can be chosen.

Non-Structured Noise Model l(n)

The non-structured interfering noise model l(n) can be created by one of the embodiments described hereinbelow.

According to an embodiment of the invention, the spectrum of the non-structured interfering noise l(n) is defined by a model having a distribution of probability of Gaussian type.

According to an embodiment of the invention, the spectrum of the non-structured interfering noise l(n) is defined by a model having a distribution of normal complex probability NC(0, 2 σ⁻¹) of zero average and of variance 2σ⁻¹, for example according to the following equation:

[Math. 50]

∀n∈{1, . . . ,N}[F _(N) l]_(n) |σ˜NC(0,2 σ⁻¹)

It can be that the background noise is non-stationary. In this case, the procedure can be via temporal segments j as was done for the interfering signal e(n). Fourier transforms S_(j), limited to a temporal segment of size N, for l(n), are defined as follows:

[Math. 51]

∀j∈{1, . . . ,J}∀n∈{1, . . . ,N _(j)}[S _(j) l]_(n)|β_(j),ζ_(j) ˜NC(0,2τ_(j) ⁻¹)

where j, J, n, N_(j), S_(j), β_(j), ζ_(j) have the definitions mentioned hereinabove,

NC(0, 2τ_(j) ⁻¹) designates a normal complex distribution of zero average, whereof the variance is equal to 2τ_(j) ⁻¹.

In each temporal window j, the non-structured interfering noise l(n) is Gaussian and has a variance 2τ_(j) ⁻¹. This is shown as

[Math. 52]

[P _(j) l]_(n)|τ_(j) ˜NC(0,2τ_(j) ⁻¹)

Total noise model b(n)

The total interfering noise model b(n) can be created by one of the embodiments described hereinbelow.

According to an embodiment of the invention, the spectrum of the interfering noise b(n) results from the models of the structured interfering noise and from the non-structured noise model. In fact, since e(n) and l(n) are independent, the density of probability of the sum b (n)=e (n)+l(n) is the convolution of the density of probability of the structured interfering noise e(n) with that of the non-structured noise l(n).

According to an embodiment of the invention, for example, in the case of a Bernoulli-Gaussian law for the structured noise, the spectrum of the total interfering noise b(n) is defined by a model having a mix of two distributions of Gaussian type each having an associated weight and different variances.

According to an embodiment of the invention, the spectrum of the total interfering noise b(n) is defined by a model having a distribution given by the barycentre between the distribution NC(0, 2τ_(j) ⁻¹ of Gaussian type weighted by the weight 1−β_(j) and a distribution NC(0, 2τ_(j) ⁻¹+τ_(j) ⁻¹)) of Gaussian type weighted by the weighting β_(j), for example according to the following equation:

[Math. 53]

∀j∈{1, . . . ,J},∀n∈{1, . . . ,N _(j)}[S _(j) b]_(n)|τ_(j),ζ_(j),β_(j) ˜(1−β_(j))NC(0,2τ_(j) ⁻¹)+β_(j) NC(0,2(ζ_(j) ⁻¹+τ_(j) ⁻¹))

where j, J, n, N_(j), S_(j), β_(j), ζ_(j), τ_(j) have the definitions mentioned hereinabove,

NC(0, 2(ζ_(j) ⁻¹+τ_(j) ⁻¹)) designates a normal complex distribution of zero average, whereof the variance is equal to 2 (ζ_(j) ⁻¹+τ_(j) ⁻¹).

As follows,

[Math. 54]

η_(j)=1/(ζ_(j) ⁻¹+τ_(j) ⁻¹)

Therefore, in this embodiment, the spectrum of the total interfering noise b(n) is defined by a model having a distribution given by the barycentre between the distribution NC(0, 2τ_(j) ⁻¹) of Gaussian type weighted by the weight 1−β_(j) and a distribution NC(0,2η_(j) ⁻¹) of Gaussian type weighted by the weight) β_(j), for example according to the following equation:

[Math. 55]

∀j∈{1, . . . ,J}∀n∈{1, . . . ,N _(j)}[S _(j) b]_(n)|τ_(j),ζ_(j),β_(j)˜(1−β_(j))NC(0,2τ_(j) ⁻¹)+β_(j) NC(0,2η_(j) ⁻¹)

where j, J, n, N_(j), S_(j), β_(j), ζ_(j), τ_(j), η_(j) have the definitions mentioned hereinabove,

NC(0,2η_(j) ⁻¹) designates a normal complex distribution of zero average, whereof the variance is equal to 2_(η) _(j) ⁻¹,

with in η_(j)<τ_(j).

τ_(j) and η_(j) are also called noise precision parameters.

The law of observation y^(R), y^(I) of the measurement signal y(t) is given by the Gaussian N according to the following hierarchical model:

$\begin{matrix} \left\{ \begin{matrix} \begin{matrix} {{{\forall{j \in {\left\{ {1,\ldots,J} \right\}{\forall{n \in {\left\{ {1,\ldots,N} \right\} d_{j}(n)}}}}}}❘\tau_{j}},} \\ {\eta_{j},{{\left. \beta_{j} \right.\sim\left( {1 - \beta} \right)\delta_{\tau_{j}^{- 1}}} + {\beta_{j}\delta_{\eta_{j}^{- 1}}}}} \end{matrix} \\ {{y^{R}❘a},{\left. d \right.\sim{N\left( {\lbrack{Ca}\rbrack^{R},{2{\sum\limits_{j = 1}^{J}{S_{j}^{\top}{{diag}\left( d_{j} \right)}S_{j}}}}} \right)}}} \\ {{y^{I}❘a},{\left. d \right.\sim{N\left( {\lbrack{Ca}\rbrack^{I},{2{\sum\limits_{j = 1}^{J}{S_{j}^{\top}{{diag}\left( d_{j} \right)}S_{j}}}}} \right)}}} \end{matrix} \right. & \left\lbrack {{Math}.56} \right\rbrack \end{matrix}$

where

the matrix C is the diagonal matrix formed by the terms exp(jΦ_(k)(n)) for n ranging from 1 to N and k ranging from 1 to K,

a is the vector of the envelopes,

diag(d_(j)) is a matrix the diagonal of which is formed by the coefficients d_(j) and whereof the other coefficients are zero,

d_(j) is a variable which follows a distribution given by the barycentre between a Dirac distribution

[Math. 57]

δ_(τ) _(j) ⁹⁻¹

positioned in the abscissa τ_(j) ⁻¹ and weighted by the weight 1−β_(j) and a Dirac distribution

[Math. 58]

δ_(η) _(j) ⁻¹

positioned in the abscissa η_(j) ⁻¹ and weighted by the weight β_(j).

According to an embodiment of the invention, the law a posteriori of d_(j) is a distribution given by the barycentre between a Dirac distribution

[Math. 59]

δ_(τ) _(j) ⁻¹

positioned in the abscissa τ_(j) ⁻¹ and weighted by the weight 1-p_(j)(n) and a Dirac distribution

[Math. 60]

δ_(η) _(j) ⁻¹

positioned in the abscissa η_(j) ⁻¹ and weighted by the weight p_(j)(n), for example according to the following equation: P [Math. 61]

∀j∈{1, . . . ,J}

[Math. 62]

∀n∈{1, . . . ,N _(j)}

[Math. 63]

d _(j)(n)|τ_(j),η_(j),β_(j) a˜(1−p _(j)(n))δ_(τ) _(j) ⁻¹ +p _(j)(n)δ_(η) _(j) ⁻¹ +p _(j)(n)δ_(η) _(j) ⁻¹

In the above,

$\begin{matrix} {{p_{j}(n)} = \frac{w_{j}(n)}{1 + {w_{j}(n)}}} & \left\lbrack {{Math}.64} \right\rbrack \end{matrix}$ $\begin{matrix} {{w_{j}(n)} = {{\frac{\beta_{j}}{1 - \beta_{j}}{\exp\left( {{- \frac{\left( {\eta_{j} - \tau_{j}} \right)\eta_{j}}{2\tau_{j}}}{❘\left\lbrack {S_{j}\left( {{Ca} - y} \right)} \right\rbrack_{n}❘}^{2}} \right)}{if}n} \in \left\{ {1,N_{j}} \right\}}} & \left\lbrack {{Math}.65} \right\rbrack \end{matrix}$ $\begin{matrix} {{w_{j}(n)} = {\frac{\beta_{j}}{1 - \beta_{j}}{\exp\left( {{- \frac{\left( {\eta_{j} - \tau_{j}} \right)\eta_{j}^{2}}{2\tau_{j}^{2}}}{❘\left\lbrack {S_{j}\left( {{Ca} - y} \right)} \right\rbrack_{n}❘}^{2}} \right)}{if}{not}}} & \left\lbrack {{Math}.66} \right\rbrack \end{matrix}$

Also, diag(d_(j)) or diag(d_(j) (n)) is a diagonal covariance matrix of the noise.

Iteration of Calculation of the Complex Envelopes a_(k)(n) of the Sinusoidal Components

According to the invention, the monitoring device 1 comprises a calculation module 12 configured, at each respective iteration carried out, to:

-   -   during a first step C1 of the respective iteration carried out         update the laws (calculate or give a value to the parameters         (average and variance) of the law) of the complex envelopes         a_(k)(n) of the sinusoidal components,     -   during a second step C21 of the respective iteration carried out         (after the first step C1 of the respective iteration carried         out) update the parameters of the noise model b(n) and/or e(n)         and/or l(n) of the measurement signal (y(t)) from the updated         complex envelopes a_(k)(n) or given the updated complex         envelopes a_(k)(n),     -   during a third step C4 of the respective iteration carried out         (after the second step C21 of the respective iteration carried         out) test whether the model converges from the preceding         iteration (t−1) to the current iteration (t), to:         -   in the negative (NO in FIG. 2 ), redo a new iteration of the             steps C1, C21, C4 described hereinabove,     -   and in the affirmative (YES in FIG. 2 ), during the fourth step         D following the third step C4, perform a calculation of the         complex envelopes a_(k)(n) from the iterations carried out, to         obtain automatic extraction of the sinusoidal components         a_(k)(n).

The updating of the noise parameters can be done independently on the j ranging from 1 to J. At each iteration t, given the laws a posteriori of the variables at the preceding iteration t−1 and the observation y, laws a posteriori of the different variables are updated by the calculation module 12, for example according to the diagram shown in FIG. 3 .

The third test step C4 enables to verify whether the laws of the different variables of the problem to be estimated at a given iteration are definitely optimal. According to an embodiment of the invention, the test of the third step C4 comprises automatic calculation of one or more indicators by the calculation module 12 from the envelopes a_(k)(n) having been calculated (for example their absolute value).

According to an embodiment of the invention, the test of the third step C4 can be performed by the calculation module 12 which compares the absolute value of the difference between the average of the law of the envelope a_(k)(n) calculated from the current iteration and the average of the law of the envelope a_(k)(n) calculated from the preceding iteration to a threshold prescribed in advance, to determine whether this difference is less than this threshold (affirmative convergence hereinabove (YES in FIG. 2 ) when this difference is less than the threshold; negative convergence hereinabove (NO in FIG. 2 ) when this difference is not less than the threshold). For example, this threshold can be prescribed between 10⁻³ and 10⁻⁶ times the norm of the average of the law of the envelope a_(k) (n) calculated from the preceding or current iteration. The test of the third step C4 can also be performed by the calculation module 12 on one or more of the described variables of the noise model, such as for example an average, a variance or other, for example of the calculated envelope a_(k) (n), to test if this variable converges from the preceding iteration to the current iteration. As long as the statistics are not stabilised (NO in FIG. 2 ), the algorithm continues to update the variables iteratively.

According to an embodiment of the invention, during the fourth step D, the calculation of the complex envelopes a_(k) (n) can be performed by the calculation module 12 from the average of the estimated law. Additional statistics can also be given as the confidence interval.

According to an embodiment of the invention, a maximal number of iterations is prescribed in advance in the calculation module 12.

Parameters of the Noise Model

According to an embodiment of the invention, these parameters can be one or more of the parameters configuring the Gaussian law of the non-structured noise (average, precision) and those configuring the heavy-tailed law of the structured noise (average, precision, form parameter etc) and this for one or more windows of the J windows.

According to an embodiment of the invention, in the case of a Bernoulli-Gaussian model these parameters can be one or more or all of the τ_(j) for j ranging from 1 to J, and/or one or more or all of the β_(j) for j ranging from 1 to J, and/or one or more or all of the η_(j) for j ranging from 1 to J, and/or one or more or all of the ζ_(j) for j ranging from 1 to J.

This embodiment of the parameters can be realised by one or more of the other embodiments described hereinbelow.

According to an embodiment of the invention, a law of probability a priori is predefined for these parameters.

According to an embodiment of the invention, conjugated laws are selected. This facilitates calculations by the calculation module 12.

According to an embodiment of the invention, in the case of a Bernoulli-Gaussian model for structured noise, a law of uniform probability over the interval ranging from 0 to 1 is allocated to one or more or all the β_(j), that is

[Math. 67]

β_(j) ˜U(0,1)∀j∈{1, . . . ,J}

According to an embodiment of the invention, a law of probability of Beta law type is allocated to one or more or all the β_(j).

By way of reminder, by convention and in general, a Beta law (noted Beta(α, H)) having a third adjustment α (strictly positive real) and a fourth adjustment H (strictly positive real) is defined by the function of density of probability f(x; α, H) according to the equation hereinbelow:

$\begin{matrix} {{f\left( {{x;\alpha},H} \right)} = \left\{ \begin{matrix} {{\frac{{x^{\alpha - 1}\left( {1 - x} \right)}^{H - 1}}{\int_{0}^{1}{{u^{\alpha - 1}\left( {1 - u} \right)}^{H - 1}{du}}}{for}x} \in \left\lbrack {0,1} \right\rbrack} \\ {0{if}{not}} \end{matrix} \right.} & \left\lbrack {{Math}.68} \right\rbrack \end{matrix}$

According to an embodiment of the invention, the law a posteriori of one or more or all the β_(j) resulting from a uniform law a priori and given the observations y(n) is a law Beta(2n₁+1, 2n₂+1) according to the following equation:

$\begin{matrix} {\forall{j \in \left\{ {1,\ldots,J} \right\}}} & \left\lbrack {{Math}.69} \right\rbrack \end{matrix}$ $\begin{matrix} {{\beta_{j}❘y},a,d,\tau_{j},{\left. \eta_{j} \right.\sim{{Beta}\left( {{{2.n_{1}} + 1},{{2.n_{2}} + 1}} \right)}}} & \left\lbrack {{Math}.70} \right\rbrack \end{matrix}$ with $\begin{matrix} {n_{1} = {\left\langle {{d_{j}(1)} = \eta_{j}^{- 1}} \right\rangle + \left\langle {{d_{j}(N)} = \eta_{j}^{- 1}} \right\rangle + {\sum\limits_{n = 2}^{N_{j} - 1}\left\langle {{d_{j}(n)} = \eta_{j}^{- 1}} \right\rangle}}} & \left\lbrack {{Math}.71} \right\rbrack \end{matrix}$ $\begin{matrix} {n_{2} = {\left\langle {{d_{j}(1)} = \tau_{j}^{- 1}} \right\rangle + \left\langle {{d_{j}(N)} = \tau_{j}^{- 1}} \right\rangle + {2{\sum\limits_{n = 2}^{N_{j} - 1}\left\langle {{d_{j}(n)} = \tau_{j}^{- 1}} \right\rangle}}}} & \left\lbrack {{Math}.72} \right\rbrack \end{matrix}$

where the function < > used hereinabove means in general <P>=1 if P is true, or if not <P>=0.

This law of probability Beta(2n₁+1,2n₂+1) is the law a posteriori which result from the law a priori (uniform law) and from the law of observations (law of y).

According to an embodiment of the invention, a law of Gamma probability is allocated to one or more or all the τ_(j), for example

[Math. 73]

τ_(j)˜Gamma(a _(τ) _(j) ,b _(τ) _(j) )

having

[Math. 74]

a _(τ) _(j)

as first adjustment and

[Math. 75]

b _(τ) _(j)

as second adjustment.

By way of reminder, by convention and in general, a Gamma law (noted Gamma(H, α)) having a first adjustment H (strictly positive real) and a second adjustment α (strictly positive real) is defined by the function of density of probability f(x; H, α) as per the equation hereinbelow:

$\begin{matrix} {{f\left( {{x;H},\alpha} \right)} = \frac{x^{H - 1}e^{{- \alpha}x}}{{\Gamma(H)}\alpha^{H}}} & \left\lbrack {{Math}.76} \right\rbrack \end{matrix}$

where Γ(H) is the gamma Euler function according to the equation hereinbelow

[Math. 77]

Γ(H)=∫₀ ^(+∞) t ^(H−1) e ^(−t) dt

According to an embodiment of the invention, a law of Gamma probability is allocated to one or more or all the η_(j), for example

[Math. 78]

η_(j)˜Gamma(a _(η) _(j) ,b _(η) _(j) )

having

[Math. 79]

a _(η) _(j) ,

as first adjustment and

[Math. 80]

b _(η) _(j) ,

as second adjustment.

The adjustments

[Math. 81]

a _(τ) _(j) ,

[Math. 82]

b _(τ) _(j) ,

[Math. 83]

a _(η) _(j) ,

[Math. 84]

b _(η) _(j) ,

of Gamma laws are each strictly positive.

According to an embodiment of the invention,

[Math. 85]

a _(τ) _(j) ,

[Math. 86]

b _(τ) _(j) ,

[Math. 87]

a _(η) _(j) ,

[Math. 88]

b _(η) _(j) .

of the Gamma laws are each prescribed in advance in the calculation module 12. Therefore,

[Math. 89]

a _(τ) _(j) ,

[Math. 90]

b _(τ) _(j) ,

[Math. 91]

a _(η) _(j) ,

[Math. 92]

b _(η) _(j) .

of the Gamma laws are each a hyper-parameter quantifying information a priori on the first parameters to be estimated.

According to an embodiment of the invention, each couple (a, b) of first adjustment and second adjustment, such as for example

[Math. 93]

(a,b)∈{(a _(η) _(j) ,b _(η) _(j) ))_(j∈{1, . . . ,J}),(a _(τ) _(j) ,b _(τ) _(j) )_(j∈{1, . . . ,J})},

is prescribed in advance in the calculation module 12 according to information a priori on the distribution of the first parameter B considered with

[Math. 94]

θ∈{η_(j∈{1, . . . ,J}),τ_(j∈{1, . . . ,J})}.

These adjustments can be prescribed in advance by the user. If the user does not prescribe these adjustments, these adjustments are fixed to very low values close or equal to zero (see the example hereinbelow) to have non-informative laws. This information can be an average of the law considered for θ and/or a variance of the law considered for θ. According to an embodiment of the invention, if the value of θ noted by θ_(e)>0 is known approximately with uncertainty equal to εθ_(e) where ε>0, then

$\begin{matrix} {{a = \frac{\theta_{e}}{\varepsilon}}{and}} & \left\lbrack {{Math}.95} \right\rbrack \end{matrix}$ $\begin{matrix} {b = \frac{1}{\varepsilon}} & \left\lbrack {{Math}.96} \right\rbrack \end{matrix}$

are prescribed in advance. If there is no information a priori on the parameter θ, alors a=b=0 are prescribed in advance (the law is non-informative). The presence of information a priori on the parameters most often accelerates convergence.

According to an embodiment of the invention, the law a posteriori of one or more or all the η_(j) resulting from its Gamma law a priori of parameters

(

[Math. 97]

a _(η) _(j) ,

[Math. 98]

b _(η) _(j) ,

and the measurement signal y(n) is a law

[Math. 99]

Gamma(ã _(η) _(j) ,{tilde over (b)} _(η) _(j) )

having

[Math. 100]

ã _(η) _(j)

as first adjustment and

[Math. 101]

{tilde over (b)} _(η) _(j)

as second adjustment according to the following equation:

[Math. 102]

∀j∈{1, . . . ,J}

[Math.103]

η_(j) |y,a—Gamma(ã _(η) _(j) ,{tilde over (b)} _(η) _(j) )

According to an embodiment of the invention, the law a posteriori of one or more or all the τ_(j) resulting from its Gamma law a priori of parameters (

[Math.104]

a _(τ) _(j) ,

[Math. 105]

b _(τ) _(j) )

and the measurement signal y(n) is a Gamma law

[Math.106]

Gamma(ã _(τ) _(j) ,{tilde over (b)} _(τ) _(j) )

having

[Math. 107]

ã _(τ) _(j)

as first adjustment and

[Math. 108]

{tilde over (b)} _(τ) _(j) ,

as second adjustment according to the following equation:

$\begin{matrix} {\forall{j \in \left\{ {1,\ldots,J} \right\}}} & \left\lbrack {{Math}.109} \right\rbrack \end{matrix}$ $\begin{matrix} {\tau_{j}{❘{y,{a \sim {Gamma}\left( {{\overset{˜}{a}}_{\tau_{j}},{\overset{˜}{b}}_{\tau_{j}}} \right){In}{the}{above}},}}} & \left\lbrack {{Math}.110} \right\rbrack \end{matrix}$ $\begin{matrix} {{\overset{˜}{a}}_{\eta_{j}} = {a_{\eta_{j}} + n_{1}}} & \left\lbrack {{Math}.111} \right\rbrack \end{matrix}$ $\begin{matrix} {{\overset{˜}{b}}_{\eta_{j}} = {b_{\eta_{j}} + {\frac{\eta_{j}}{2}\ {\sum\limits_{\underset{{d_{j}(n)} = \eta_{j}^{- 1}}{n}}\ {❘\left\lbrack {S_{j}\left( {{Ca} - y} \right)} \right\rbrack_{n}❘}^{2}}}}} & \left\lbrack {{Math}.112} \right\rbrack \end{matrix}$ $\begin{matrix} {{\overset{˜}{a}}_{\tau_{j}} = {a_{\tau_{j}} + n_{2}}} & \left\lbrack {{Math}.113} \right\rbrack \end{matrix}$ $\begin{matrix} {{\overset{˜}{b}}_{\tau_{j}} = {b_{\tau_{j}} + {\frac{\tau_{j}}{2}\ {\sum\limits_{\underset{{d_{j}(n)} = \tau_{j}^{- 1}}{n}}\ {❘{\left\lbrack {S_{j}\left( {{Ca} - y} \right)} \right\rbrack_{n}❘^{2}}}}}}} & \left\lbrack {{Math}.114} \right\rbrack \end{matrix}$

Regularisation Parameters γ_(k)

According to an embodiment of the invention, the calculation module 12 is configured, during a fifth step C22 of the respective iteration carried out (following the second step C21 of the respective iteration carried out and prior to the third step C4) to update second regularisation parameters γ_(k) of the model from the updated complex envelopes a_(k)(n) or given the updated complex envelopes a_(k)(n). These second regularisation parameters γ_(k) are added restrictions imposed to favour regularity of the complex envelopes a_(k)(n). These second regularisation parameters γ_(k) characterise suppression of singularities in the components a_(k)(n).

This embodiment of the second regularisation parameters γ_(k) can be realised by one of the other embodiments described hereinbelow.

Updating of the second parameters can be done independently on the k ranging from 1 to K.

According to an embodiment of the invention, the density p of probability (a priori) of the complex envelopes a_(k)(n) in the presence of these second regularisation parameters γ_(k) is the following, for k ranging from 1 to K:

$\begin{matrix} {{p\left( {a^{R},{a^{I}{❘{\gamma_{1},\ldots,\gamma_{K}}}}} \right)} = {{\prod\limits_{k = 1}^{K}{{p\left( {a_{k}^{R}{❘\gamma_{k}}} \right)}{p\left( {a_{k}^{I}{❘\gamma_{k}}} \right)}}} = {C_{a}{\prod\limits_{k = 1}^{K}{\gamma_{k}^{N}{\exp\left( {{- \gamma_{k}}{{A_{p}a_{k}}}^{2}} \right)}}}}}} & \left\lbrack {{Math}.115} \right\rbrack \end{matrix}$

where C_(a) is a constant independent from the γ_(k) and from the complex envelopes a_(k)(n).

According to an embodiment of the invention, a law of Gamma probability is allocated to one or more or all the second regularisation parameters γ_(k), for example

[Math. 116]

γ_(k)˜Gamma(a _(β) _(k) ,b _(γ) _(k) )

having

[Math. 117]

a _(γ) _(k)

as first adjustment and

[Math. 118]

b _(γ) _(k)

as second adjustment. According to an embodiment of the invention, the couple (

[Math. 119]

a _(γ) _(k) ,

[Math. 120]

b _(γ) _(k) )

is prescribed in advance in the calculation module 12 according to information a priori on the distribution of the first parameter γ_(k). This information can be an average of the law considered for γ_(k) and/or a variance of the law considered for γ_(k). These adjustments can be prescribed in advance by the user. If the user does not prescribe these adjustments, these adjustments are fixed to very low values close or equal to zero (see the example hereinbelow) to have non-informative laws. According to an embodiment of the invention, if the value of γ_(k) noted by θ_(e)>0 is known approximately with uncertainty equal to εθ_(e) where ε>0, then

$\begin{matrix} {{a_{\gamma_{k}} = \frac{\theta_{e}}{\varepsilon}}{and}} & \left\lbrack {{Math}.121} \right\rbrack \end{matrix}$ $\begin{matrix} {b_{\gamma_{k}} = \frac{1}{\varepsilon}} & \left\lbrack {{Math}.122} \right\rbrack \end{matrix}$

are prescribed in advance. If there is no information a priori on the parameter θ, then

[Math. 123]

a _(γ) _(k) =b _(γ) _(k) =0

are prescribed in advance (the law is non-informative). The presence of information a priori on the parameters most often accelerates convergence.

According to an embodiment of the invention, a law of probability a posteriori of γ_(k) (resulting from the Gamma law a priori of adjustment parameters (

[Math. 124]

a _(γ) _(k)

[Math. 125]

b _(γ) _(k) ),

and from the law a priori of envelopes) is also a law Gamma

[Math. 126]

Gamma(ã _(γ) _(k) ,{tilde over (b)} _(γ) _(k) )

having

[Math. 127]

ã _(γ) _(k)

as first adjustment and

[Math. 128]

{tilde over (b)} _(γ) _(k)

as second adjustment, for example according to the following equation:

[Math. 129]

∀k ∈ {1, …, K} $\begin{matrix} {\gamma_{k}{❘{{a_{k} \sim {Gamma}\left( {{\overset{˜}{a}}_{\gamma_{k}},{\overset{˜}{b}}_{\gamma_{k}}} \right){In}{the}{above}},}}} & \left\lbrack {{Math}.130} \right\rbrack \end{matrix}$ $\begin{matrix} {{\overset{˜}{a}}_{\gamma_{k}} = {a_{\gamma_{k}} + N}} & \left\lbrack {{Math}.131} \right\rbrack \end{matrix}$ $\begin{matrix} {{\overset{˜}{b}}_{\gamma_{k}} = {b_{\gamma_{k}} + {\frac{1}{2}{{A_{p}a_{k}}}^{2}}}} & \left\lbrack {{Math}.132} \right\rbrack \end{matrix}$

According to an embodiment of the invention, a law a posteriori of the sinusoidal components of interest a_(k)(n)exp(jΦ_(k)(n)) is calculated during the first step C1 by the calculation module 12 from of the law of observation of the noise b (n) and from the law of observation of the measurement signal y(t). According to an embodiment of the invention, the law a posteriori of the real part a^(R) of the sinusoidal components of interest a_(k)(n)exp(jΦ_(k)(n)) is defined by a model having a distribution of probability of Gaussian type N(m_(x) ^(R),Σ_(x)) of average m_(x) ^(R) (real part of m_(x)) and of variance Σ_(x), according to the following equation:

[Math. 133]

a ^(r) |y,d,γ ₁, . . . ,γ_(k) ˜n(m _(x) ^(R),Σ_(x))

According to an embodiment of the invention, the law a posteriori of the imaginary part a^(I) of the sinusoidal components of interest a_(k)(n)exp(jΦ_(k)(n)) is defined during the first step C1 by a model having a distribution of probability of Gaussian type N(m_(x) ^(I),Σ_(x)) of average m_(x) ^(I) (imaginary part of m_(x)) and of variance Σ_(x), for example according to the following equation:

[Math. 134]

a ^(I) |y,d,γ ₁, . . . ,γ_(K) ˜N(m _(x) ^(I),Σ_(x))

The real part a^(R) of the sinusoidal components of interest a_(k)(n)exp(jΦ_(k)(n)) and the imaginary part a^(I) of the sinusoidal components of interest a_(k)(n)exp(jΦP_(k)(n)) can be updated independently during the first step C1. In the case of a very large number of components to be estimated, this first step C1 can be done in parallel because of multi-processor architecture of the module 12.

In the above,

$\begin{matrix} {m_{x} = {\sum_{x}{C^{\top}{\sum\limits_{j = 1}^{J}{S_{j}^{\top}{{diag}\left( d_{j} \right)}^{- 1}S_{j}y}}}}} & \left\lbrack {{Math}.135} \right\rbrack \end{matrix}$ $\begin{matrix} {\sum_{x}{= \left( {{C^{\top}{\sum\limits_{j = 1}^{J}{S_{j}^{\top}{{diag}\left( d_{j} \right)}^{- 1}S_{j}C}}} + {A^{T}{{diag}\left( \gamma_{k} \right)}A}} \right)^{- 1}}} & \left\lbrack {{Math}.136} \right\rbrack \end{matrix}$

diag(d_(j)) is a matrix whereof the diagonal is formed by the coefficients d_(j) and whereof the other coefficients are zero,

diag(γ_(k)) is a matrix whereof the diagonal is formed by the coefficients γ_(k) and whereof the other coefficients are zero,

A_(p) is the matrix of the discrete differences ∇^(p) of order p of diag(γ_(k)),

A is the diagonal block matrix formed by the K matrices A_(p) and is the matrix containing K blocks of A_(p).

the matrix C is the diagonal matrix formed by the terms exp(jΦ_(k)(n)) for n ranging from 1 to N.

Auxiliary Decoupling Variable v of the Envelopes

According to an embodiment of the invention, during the sixth step C3 (following the second step C21 of the respective iteration carried out and/or at the fifth step C22 and prior to the third step C4) an auxiliary decoupling variable v of the complex envelopes a_(k)(n) of the components between them in the model is calculated by the calculation module 12, from the updated complex envelopes a_(k)(n), from the first parameters (which can be one or more or all of the τ_(j) for j ranging from 1 to J, and/or one or more or all of the β_(j) for j ranging from 1 to J, and/or one or more or all of the η_(j) for j ranging from 1 to J, and/or one or more or all of the ζ_(j) for j ranging from 1 to J) and from the second regularisation parameters γ_(k). A known problem of optimisation algorithms is the inversion of a full matrix of large dimension, most often badly conditioned. This embodiment resolves this problem by decoupling the components by addition of the auxiliary variable v. The auxiliary variable v takes into account information on the correlations between the components a_(k)(n). Due to the auxiliary variable v, the components are estimated independently. Coupling occurs implicitly via the auxiliary variable.

This embodiment resolves the large-dimension digital problem, linked to the above first problem. This embodiment is adapted all the more since a large number K of components is preferred. In fact, the known Vold-Kalman method suffers from the problem of large-dimension matrix inversion above all when the matrix in question is very badly conditioned which is the case for very large values of r_(k). This inversion could occur in the known Vold-Kalman method due to an iterative algorithm of conjugated gradient type. This embodiment is not contradicted by this digital problem, given that the different components are decoupled in the algorithm and therefore can be estimated disjointly. The coupling between the different components is considered implicitly by the auxiliary variable v added in the updating of each component.

This embodiment of the invention on the auxiliary decoupling variable v can be realised by one of the other embodiments described hereinbelow.

According to an embodiment of the invention, the law of the auxiliary decoupling variable v is defined by a model having a distribution of probability of complex circularly symmetrical normal law type CNC(m_(v), Σ_(v)) having the average m_(v) and the variance Σ_(v). For v there is a conditional law

[Math. 137]

v|a,d˜CNC(m _(v),Σ_(v))

The complex circularly symmetrical normal law CNC(m_(v), Σ_(v)) has the average

$\begin{matrix} \begin{pmatrix} m_{v}^{R} \\ m_{v}^{I} \end{pmatrix} & \left\lbrack {{Math}.138} \right\rbrack \end{matrix}$

and covariance matrix

$\begin{matrix} {\frac{1}{2}\begin{pmatrix} {\sum_{v}^{R};} & {- \sum_{v}^{I}} \\ {\sum_{v}^{I};} & \sum_{v}^{R} \end{pmatrix}} & \left\lbrack {{Math}.139} \right\rbrack \end{matrix}$

and is defined by

$\begin{matrix} {v \sim {{CNC}\left( {m_{v},\sum_{v}} \right)}\text{<=>}\begin{pmatrix} v^{R} \\ v^{I} \end{pmatrix} \sim {N\left( {\begin{pmatrix} m_{v}^{R} \\ m_{v}^{I} \end{pmatrix},{\frac{1}{2}\begin{pmatrix} {\sum_{v}^{R};} & {- \sum_{v}^{I}} \\ {\sum_{v}^{I};} & \sum_{v}^{R} \end{pmatrix}}} \right)}} & \left\lbrack {{Math}.140} \right\rbrack \end{matrix}$

In the above,

$\begin{matrix} {m_{v} = {\sum_{v}a}} & \left\lbrack {{Math}.141} \right\rbrack \end{matrix}$ $\begin{matrix} {\sum_{v}{= {{\frac{1}{\mu}{Id}} - {C^{\top}{\sum\limits_{j = 1}^{J}{S_{j}^{\top}{{diag}\left( d_{j} \right)}^{- 1}S_{j}C}}}}}} & \left\lbrack {{Math}.142} \right\rbrack \end{matrix}$

with μ being a real prescribed non-zero verifying

$\begin{matrix} {\mu < \frac{1}{K{\max\left( {d_{j}(n)} \right)}}} & \left\lbrack {{Math}.143} \right\rbrack \end{matrix}$

This enables to eliminate problems of inversion of large-dimension matrices in particular for large values of N and K and allows an increase in data, which allows decoupling between the different components a_(k) (n)exp(jΦ_(k) (n)).

Calculation of the Sinusoidal Components a_(k) of Interest During the First Step C1

According to an embodiment of the invention, the law a posteriori of the real part a_(k) ^(R) of the sinusoidal components of interest a_(k) (n)exp(jΦ_(k) (n)) is defined by a model having a distribution of probability of normal law type having the average m_(k) ^(R) and the variance Σ_(k), for example according to the following equation:

[Math. 144]

∀k∈{1, . . . ,K}a _(k) ^(R) |y,d,γ _(k) ,v _(k) ˜N(m _(k) ^(R),Σ_(k))

According to an embodiment of the invention, the law a posteriori of the imaginary part a_(k) ^(I) of the sinusoidal components of interest a_(k)(n)exp(jΦ_(k)(n)) is defined by a model having a distribution of probability of normal law type having the average m_(k) ^(I) and the variance Σ_(k), for example according to the following equation:

[Math. 145]

∀k∈{1, . . . ,K}a _(k) ^(I) |y,d,γ _(k) ,v _(k) ˜N(m _(k) ^(I),Σ_(k))

The definitions hereinabove of a_(k) ^(R) and of a_(k) ^(I) are laws a posteriori of the components.

In the above,

$\begin{matrix} {m_{k} = {\sum_{k}\left( {{C_{k}^{\top}{\sum\limits_{j = 1}^{J}{S_{j}^{\top}{{diag}\left( d_{j} \right)}^{- 1}S_{j}y}}} + v_{k}} \right)}} & \left\lbrack {{Math}.146} \right\rbrack \end{matrix}$ $\begin{matrix} {\sum_{k}{= \left( {{\frac{1}{\mu}{Id}} + {\gamma_{k}A_{p}^{\top}A_{p}}} \right)^{- 1}}} & \left\lbrack {{Math}.147} \right\rbrack \end{matrix}$

Therefore, a_(k) ^(R) and a_(k) ^(I) can be estimated independently. The coupling is considered only via the auxiliary variable v occurring in the average m_(k). It is interesting to note that even if there is a need to make the inversion here, this operation is less expensive. On the one hand, the dimension of the matrix Σ_(k) is K times smaller than the matrix Σ_(x). On the other hand, if the choice is made for conditions of circulating edges (or zero edges by making padding in the zero edges) in the matrix of differences A_(p), the matrix Σ_(k) is then circulating and therefore diagonalisable in the Fourier field.

According to an embodiment of the invention,

$\begin{matrix} {\sum_{k}{= {{F_{N}^{\top}\left( {{\frac{1}{\mu}{Id}} + {\gamma_{k}Q_{k}^{\top}Q_{k}}} \right)}^{- 1}F_{N}}}} & \left\lbrack {{Math}.148} \right\rbrack \end{matrix}$

where Q_(k) is the diagonal matrix containing the Fourier coefficients of the filter associated with the discrete matrix of difference A_(p) which is equal to the first column of A_(p). Updating of each m_(k) can also be done by shifting to the Fourier field. Also, there is no need to store Σ_(k), or its inverse, but simply the diagonal values of the matrix

$\begin{matrix} {{\frac{1}{\mu}{Id}} + {\gamma_{k}Q_{k}^{\top}Q_{k}}} & \left\lbrack {{Math}.149} \right\rbrack \end{matrix}$

According to an embodiment of the invention, the calculation module 12 is configured to perform a seventh decision step E following the fourth step D. According to an embodiment of the invention, the decision during the seventh step E is taken automatically by the calculation module 12 from analyses carried out automatically by the calculation module 12 on the extracted components a_(k)(n) having been obtained at the fourth step D, and/or from test bench vibratory characterisation (for example identification of eigen modes) carried out by the calculation module 12 on the extracted components a_(k)(n) having been obtained at the fourth step D of the rotating machine 100. According to an embodiment of the invention, the seventh step E is performed by the calculation module 12 on a monitoring decision and/or on the triggering of an alarm AL or alert AL of defects of the rotating machine 100 (for example: imbalance, meshing, misalignment, mass eccentricity, worn bearings, twisted shafts, misaligned shafts, parallel alignment defect of a rotor shaft of the rotating machine, broken shaft, premature wear of bearings, joints, shafts and couplings of the kinematic chain, resonance). According to an embodiment of the invention, the monitoring device 1 is embedded on the aircraft. According to an embodiment of the invention, the calculation module 12 automatically controls communication of the alert, of the alarm AL or of the decision it has taken at the seventh step E on the physical outlet SP, for example on the display screen, and/or at a fixed ground station SP.

Updating of the law in the proposed device can be done in two ways.

On the one hand, it can be approximated by samples. This concerns the Monte Carlo Markov Chain methods (MCMC). Therefore, updating a law means taking a sample of this law. It is noted that all laws are simple to sample. For the parameters of the law of observation and of the regularisation parameter, this is evident. For the law a posteriori of the envelopes, the generation of random variables can be done by shifting to the Fourier field. And given the particular form of Σ_(v), the generation of random auxiliary variables can also be done directly based on the fact that CC^(T)=K. Id and S_(j)S_(j) ^(T)=Id. After enough iterations, the random variables generated by the algorithm will follow the preferred laws a posteriori; therefore the statistics (average, variance etc) can be approximated by using empirical estimators with these samples.

On the other hand, the law a posteriori can be approximated by another simpler one, for example by a Bayesian-Variational (ABV). approximation approach. In this way the statistics of laws of interest can be approximated by those of approximating laws after enough iterations.

An example of application of the invention on a vibratory signal of a helicopter is given hereinbelow in FIGS. 4 a to 4 f . The aim is to monitor the component N₂ of the vibratory signal y(t) of the high-pressure shaft of a turbomachine 100, for example of the shaft of the free turbine in the helicopter engine.

FIG. 4 a shows the spectrogram of the measured signal y(t) by indicating, for each instant in abscissae and each frequency on the frequency band [0,700] Hz in ordinates, a pixel whereof the greyscale corresponds according to the right-hand scale of values ECH to a spectral amplitude for this frequency (and similarly for the spectrograms of FIGS. 4 c, 4 d, 4 e and 4 f ). FIG. 4 b shows in ordinates the amplitude of the envelope a_(k)(n) having been calculated according to the invention as a function of time in abscissae (noise is considered non-stationary and non-Gaussian and the regularisation parameter is estimated automatically). FIG. 4 b shows a significant rise in amplitude from the instant 5 s, which is the sign of a passage through a resonance. The evaluation of the maximal amplitude during this resonance is highly important practically for vibratory control of the aircraft in transient operation. FIG. 4 c shows the spectrogram of the residual signal (the measured signal y(t) from which the component N₂ is subtracted). FIG. 4 c shows that the component N₂ of interest has been correctly estimated. FIGS. 4 d, 4 e and 4 f show the spectrograms of the residual signal by using the known method Vold-Kalman with different regularisation parameters adjusted manually (these are values currently used in literature). FIG. 4 d shows an underestimation of the amplitude of the component N₂ of interest when passing through resonance for a large value of r_(k) (the estimated amplitude of the component N₂ of interest is less important than the real one, the estimated envelope is very smoothed excessively, which very aggressively penalises the great variations of the envelope). FIGS. 4 e and 4 f show that there is overestimation of the amplitude of the component N₂ of interest, for a small value of r_(k) (the estimated amplitude of the component N₂ of interest is more important than the real amplitude, the estimated envelope is very affected by noise). In fact, the component N₂ of interest having been estimated in FIGS. 4 e and 4 f is badly filtered and contains much background noise. Also, overestimation of the envelope can result from a badly conditioned matrix and therefore from digital problems (see FIGS. 1 and 4 f for example).

According to another embodiment, β_(j) is equal to zero. In this case the contribution of the narrowband noise e(n) is omitted. The principal advantage is the reduction in the number of parameters to be estimated (β_(j) and ζ_(j)) as well as the calculation cost. Admittedly, this can be a good approximation if no parasite harmonic of the residual signal interferes directly with the signals of interest. This goes back to y|a, τ˜NC (Ca,2τ⁻¹)

In this particular case, the advantage relative to the known Vold-Kalman method is the automatic adjusting of the regularisation parameter γ_(k). In fact, the Vold-Kalman parameter r_(k) is linked to the parameters of the invention by r_(k)=γ_(k)/τ.

A disadvantage of this approximation could manifest, for example, in the presence of interception between the harmonics of interest and the residual harmonics. In fact, estimation artefacts can occur at instants of interferences.

However, if the instants of interferences are known approximately (for example with inspection of the spectrogram), it can be assumed that the noise is Gaussian except for close to coupling, that is at intervals j ∈{1, . . . , J} containing the coupling instants. This means imposing β_(j)=0 in the other intervals not containing the coupling moments.

Of course, the embodiments, characteristics, possibilities and examples described hereinabove can be combined with each other or be selected independently of each other.

If a parsimonious model other than the Bernoulli-Gaussian model is used, the same forms of laws are obtained. In this case, the parameters of the noise models are going to change. The law of the envelopes, of the regularisation parameters and the auxiliary variable will have the same form. The chosen noise model will act only on the values of d_(j). 

1. A method for monitoring at least one rotating machine of an aircraft, the method comprising the following steps: acquiring at least one measurement signal of the rotating machine is by an acquisition module, estimating instantaneous frequencies of sinusoidal components of the at least one measurement signal by an estimation module, executing by a calculation module several successive iterations, each of which comprising: updating complex envelopes of the sinusoidal components, updating parameters of a noise model of the at least one measurement signal from the complex envelopes having been updated, performing a convergence test as to whether the noise model converges from the preceding iteration to the current iteration, to: in the negative, redo a new iteration, in the affirmative, perform a calculation of the complex envelopes from the iterations carried out, updating at each iteration regularisation parameters of the noise model by the calculation module from the complex envelopes having been updated, and, updating, at each iteration, an auxiliary decoupling variable of the complex envelopes of the components between them in the noise model by the calculation module from the complex envelopes having been updated and from the parameters of the noise model having been updated.
 2. The method according to claim 1, wherein the noise model of the at least one measurement signal comprises a structured interfering noise, whereof a distribution of probability a priori of spectrum has a heavy-tailed law calculated by the calculation module at each iteration.
 3. The method according to claim 1, wherein the noise model of the at least one measurement signal comprises a non-structured interfering noise, whereof a distribution of probability a priori is given by a third distribution of Gaussian type calculated by the calculation module at each iteration.
 4. The method according to claim 1, wherein the noise model of the at least one measurement signal resulting from a structured interfering noise model and from a non-structured interfering noise model is a convolution of a Gaussian law and of a heavy-tailed law whose parameters are calculated by the calculation module at each iteration.
 5. The method according to claim 1, wherein a conjugated law of probability a priori is allocated to hyperparameters of laws a priori of the parameters of the noise model.
 6. The method according to claim 1, wherein a Gamma law of probability a priori is allocated to at least one of the regularisation parameters and results from a Gamma law a posteriori whose parameters are calculated by the calculation module at each iteration.
 7. The method according to claim 1, wherein the auxiliary decoupling variable is defined by a model having a distribution of probability of a type of circularly symmetrical normal complex law, which is calculated by the calculation module at each iteration.
 8. The method according to claim 1, wherein the convergence test comprises testing by the calculation module at each current iteration that the difference between an absolute value of the complex envelope calculated from the current iteration and the absolute value of the complex envelope calculated from the preceding iteration is less than a prescribed threshold.
 9. The method according to claim 1, comprising analysing by the calculation module the complex envelopes having been calculated, and triggering a communication of an alarm by the calculation module on at least one physical outlet as a function of the complex envelopes having been analysed.
 10. A monitoring device of at least one rotating machine of an aircraft, the device comprising a module for acquiring at least one measurement signal of the rotating machine, wherein the device comprises a module for estimating instantaneous frequencies of sinusoidal components of the at least one measurement signal, a calculation module configured, at each of several successive iterations carried out, to: update complex envelopes of the sinusoidal components, update parameters of a noise model of the at least one measurement signal from the complex envelopes having been updated, test whether the noise model converges from the preceding iteration to the current iteration, to: in the negative, redo a new iteration, in the affirmative, perform a calculation of the complex envelopes from the iterations carried out, the calculation module being configured to, at each of the successive iterations carried out: update regularisation parameters of the noise model from the complex envelopes having been updated, update an auxiliary decoupling variable of the complex envelopes of the components between them in the noise model, from the complex envelopes having been updated and from the parameters of the noise model having been updated.
 11. The monitoring device according to claim 10, wherein the calculation module is configured to analyse the complex envelopes having been calculated, and is able to trigger a communication of an alarm on at least one physical outlet of the device as a function of the complex envelopes having been analysed.
 12. The monitoring device according to claim 10, wherein the auxiliary decoupling variable is defined by a model having a distribution of probability of a type of circularly symmetrical normal complex law, which is calculated by the calculation module at each iteration.
 13. A computer program comprising code instructions for executing the monitoring method according to claim 1, when it is executed on a computer.
 14. An aircraft comprising at least one rotating machine and a monitoring device of the at least one rotating machine according to claim
 10. 